{"cells":[{"cell_type":"markdown","metadata":{"id":"22B733A3BFE049CCA72BE2F628394282","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":" # <center>  Lecture 6 : Approximating the Posterior </center>  \n \n##  <center>  Instructor: Dr. Hu Chuan-Peng   </center>  \n"},{"cell_type":"markdown","metadata":{"id":"E2A68D14F12640178A433C789F678B13","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"**Intro**  \n\n我们回到某选举的例子，在Lecture 3中，选举结果只受到一个单一支持率$“\\pi”$的影响，但这是一个简化的模型。  \n\n如果我们想让这个模型更接近现实，那么更多的影响因素就会被包含在内，例如：  \n- 公众形象  \n- 军事威慑  \n- 执政能力  \n\n和支持率$pi$一样，这些影响因素都可以用参数来表示 $\\theta = (\\theta_1, \\theta_2, \\ldots, \\theta_k)$  \n\n\n![Image Name](https://cdn.kesci.com/upload/s2ldkt43rd.png?imageView2/0/w/640/h/640)  \n"},{"cell_type":"markdown","metadata":{"id":"8ED2FBC4EC7848C88EC1C96E1005E4D4","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"回想我们之前关于贝叶斯的内容：  \n\n- 单个事件的贝叶斯更新  \n- 单个变量，$\\theta$的更新  \n\n但现实生活中我们要解决的问题呢？  \n"},{"cell_type":"markdown","metadata":{"id":"9A1980E50F464841851D6638076B7884","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"现实生活中的问题往往设计多个因素(参数)，而当代入更多的参数$\\theta_k$时：  \n\n$$  \nf(\\theta_k | y) = \\frac{f(\\theta_k)L(\\theta_k | y)}{f(y)}  \n$$  \n\n* 此时的$f(\\theta_k)$与$L(\\theta_k|y)$都是多元函数，即他们返回的结果是关于多个参数而不仅是一个参数。  \n\n问题在于，但此时参数$\\theta$的数量增多，计算$f(y)$的计算难度就会变大：  \n\n$$  \nf(y) = \\int_{\\theta_1}\\int_{\\theta_2} \\cdots \\int_{\\theta_k} f(\\theta)L(\\theta | y) d\\theta_k \\cdots d\\theta_2 d\\theta_1  \n$$  \n\n* 对于一个参数$\\theta$来说，$f(y) = f(\\theta)$ * $L(\\theta|y)$是$\\theta$的所有取值下所以似然值的积分(可以理解为求和) $\\int_{\\theta} f(\\theta)L(\\theta | y) d\\theta$  \n* 当参数变多，需要对所有的参数组合$\\int_{\\theta_1}\\int_{\\theta_2} \\cdots \\int_{\\theta_k} f(\\theta)L(\\theta | y) d\\theta_k \\cdots d\\theta_2 d\\theta_1$进行积分，计算量会大大增加。  \n* 有些似然函数可能难以求解。"},{"cell_type":"markdown","metadata":{"id":"869A2261F904486B99C4020DB4D3186D","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"问题：是否可以通过一些**近似**的方法来近似后验分布？  \n\n我们将介绍两种近似方法：  \n1. 网格近似  \n2. MCMC"},{"cell_type":"markdown","metadata":{"id":"EFE827F84E5944F59E01A7C7CFD2AD62","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"## Grid approximation 网格近似  \n![Image Name](https://www.bayesrulesbook.com/bookdown_files/figure-html/unnamed-chunk-177-1.png)  \n\n\n> * 想象有一张图像(图三)，但你不能完整地看到它。不过你可以从左到右每次取一个这个图像中的一个小格来观察它。  \n> * 只要格子越细，最后组合在一起就会与完整的图片更近似  \n\n在网格近似中，后验分布$f(\\theta | y)$其实就是完整的图片。我们可以选择有限个$\\theta$，并观察对应的$f(\\theta | y)$，以此来近似完整的后验分布。  \n"},{"cell_type":"markdown","metadata":{"id":"0CFE4841C41D4DDA82B1CD7047836DAB","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"（以一个参数为例）网格近似可以分为以下四个步骤：  \n1. 选定一系列离散的$\\theta$值  \n2. 计算每个$\\theta$值对应的先验分布$f(\\theta)$和似然函数$L(\\theta|y)$  \n3. 对于所有的$\\theta$值，计算$f(\\theta)$与$L(\\theta|y)$二者的乘积并相加，再进行归一化  \n4. 归一化后，根据$\\theta$值的后验概率分布，随机抽取$N$个$\\theta$值  \n\n我们以Beta-Binomial的例子来演示这四个步骤"},{"cell_type":"markdown","metadata":{"id":"CDEF576C027248DCAF7152DC6CC3C7B2","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"### A Beta-Binomial example  \n\n现在假设先验$\\pi$ 服从  \n\n$$  \n\\pi    \\sim \\text{Beta}(2, 2)  \n$$  \n\n似然函数为：  \n$$  \nY|\\pi  \\sim \\text{Bin}(10, \\pi)  \n$$  \n\n> * 这是我们在lec3熟悉的beta-binomial模型，假设$\\pi$反映的是支持率，似然函数反映的则是在某个支持率下，总投票数为10时，支持票数的分布概率情况。  \n> $Y | \\pi  \\sim \\text{Bin}(n, \\pi)$  \n> $\\pi \\sim \\text{Beta}(\\alpha, \\beta)$  \n> * 在观察到 n 次时间中有 Y = y次目标事件后，$\\pi$的后验分布可以用Beta模型来描述，反映了先验（通过α和β）和数据（通过y和n）的影响：  \n> $\\pi | (Y = y) \\sim \\text{Beta}(\\alpha + y, \\beta + n - y) $  \n"},{"cell_type":"markdown","metadata":{"id":"A196A9EA235149CD9F73EA110FA1D0F2","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"假设支持票数为9，我们通过共轭先验的公式，可以直接计算出$Beta$后验分布的两个参数：  \n- $Y + \\alpha = 9 + 2\\;\\;, \\;\\;\\;\\;\\;\\;\\;\\; n - Y + \\beta = 10 - 9 + 2$  \n- $\\pi | (Y = 9) \\sim \\text{Beta}(11, 3)$"},{"cell_type":"markdown","metadata":{"id":"525416FB18444AD0B78BEFB872BC1C60","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"**代码示例：网格近似估计后验分布**  \n\n现在我们暂时忘记后验的简便计算方法，使用网格近似来对后验分布进行估计。  \n\n**step1：** 首先，我们将从0~1内的连续变量$\\pi$中取出6个值：$\\pi \\in \\{0, 0.2, 0.4, 0.6, 0.8, 1\\}$"},{"cell_type":"code","metadata":{"collapsed":false,"id":"790A302DD36441088FCFF87ADDB3201A","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"skip"},"tags":[],"trusted":true},"source":"import numpy as np\nimport scipy.stats as st\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport pandas as pd","outputs":[],"execution_count":2},{"cell_type":"code","metadata":{"collapsed":false,"id":"7CCD42149EFB403DA8B3863DED43A267","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"source":"pi_grid = np.linspace(0,1,6)\nprint(\"从0~1内的连续变量π中取出6个值：\", pi_grid)","outputs":[{"name":"stdout","output_type":"stream","text":["从0~1内的连续变量π中取出6个值： [0.  0.2 0.4 0.6 0.8 1. ]\n"]}],"execution_count":5},{"cell_type":"markdown","metadata":{"id":"047D875ED8D74AF884831F0455B60E48","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"**step2、3：**  \n    在每一个$\\pi$下，计算先验分布 $Beta(2, 2)$，与似然函数$Bin(10, \\pi)$ (Y=9)的乘积，计算总和并进行归一化"},{"cell_type":"code","metadata":{"collapsed":false,"id":"2ED3397AD2EF4272A2E0D45A7019C80A","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"source":"prior = st.beta.pdf(pi_grid, 2, 2)\nlikelihood = st.binom.pmf(9, 10, pi_grid)\nposterior = prior * likelihood / np.sum(prior * likelihood)\n\nplt.stem(pi_grid, posterior,    # 使用 plt.stem() 绘制垂直柱状图，表示在 pi_grid 上的 posterior 分布\n         linefmt='black',\n         bottom=-1)  \n\nplt.ylim(0, 1)                  # 设置 y 轴范围在 0 到 1 之间\nplt.xlabel('pi_grid')           # 设置 y 轴范围在 0 到 1 之间\nplt.ylabel('Posterior')         # 设置 y 轴范围在 0 到 1 之间\nsns.despine()                   # 使用 sns.despine() 函数去除图形的上方和右侧的边框","outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/5CF08C438D7049D19178F0AE1DF86F0E/s2m1lqfrug.png\">"],"text/plain":["<Figure size 640x480 with 1 Axes>"]},"metadata":{},"output_type":"display_data"}],"execution_count":7},{"cell_type":"markdown","metadata":{"id":"50717EBAAE85435CA94CF4D282BB4F93","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"**step4：**  \n    得到后验结果($\\pi$)的分布后，从这个分布中抽样10000次  \n\n*由于$\\pi$的取值只在 $\\{0, 0.2, 0.4, 0.6, 0.8, 1\\}$之中，从上图中可以看到，只有0.6与0.8的概率值相对较高"},{"cell_type":"code","metadata":{"collapsed":false,"id":"1788106586FB4C9486C45F9D83FE0F32","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"source":"np.random.seed(84735)\n\nposterior_sample = np.random.choice(pi_grid,        # pi_grid 为抽样的候选集合\n                                    size=10000,     # 从 posterior 分布中随机抽取 10000 个样本\n                                    p=posterior,    # p=posterior 为抽样的概率分布\n                                    replace=True)   # replace=True 表示有放回抽样\n\nposterior_sample = pd.DataFrame({\"pi_sample\": posterior_sample})    # 将抽取的样本存储在 DataFrame 中，列名为 \"pi_sample\"\n\nposterior_sample.value_counts(normalize=True).reset_index()        # 对 posterior_sample 中的样本进行计数，并使用 normalize=True 将计数转换为相对频率\n                                                                   # 然后通过 reset_index() 将结果重新设置为 DataFrame","outputs":[{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th>pi_sample</th>\n","      <th>0</th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>0</th>\n","      <td>0.8</td>\n","      <td>0.8060</td>\n","    </tr>\n","    <tr>\n","      <th>1</th>\n","      <td>0.6</td>\n","      <td>0.1858</td>\n","    </tr>\n","    <tr>\n","      <th>2</th>\n","      <td>0.4</td>\n","      <td>0.0082</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>"],"text/plain":["   pi_sample       0\n","0        0.8  0.8060\n","1        0.6  0.1858\n","2        0.4  0.0082"]},"execution_count":8,"metadata":{},"output_type":"execute_result"}],"execution_count":8},{"cell_type":"markdown","metadata":{"id":"A928931268E3426C94471E5A4853CF54","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"**抽样结果图示**  \n\n可以看到，在抽样的结果中，$\\pi$大多都是0.6或0.8，少部分是0.4  \n\n我们对比一下抽样得到的后验分布图 和 实际的后验分布$Beta(11, 3)$"},{"cell_type":"code","metadata":{"collapsed":false,"id":"E095D6CB259B475D9AE58B754C0A79C3","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"source":"x_beta = np.linspace(0, 1, 10000)       # 生成10000 个点，范围在 [0, 1] 之间\n\ny_beta = st.beta.pdf(x_beta, 11, 3)     # 生成 Beta(11,3)\n\nplt.plot(x_beta, y_beta, color='black') # 绘制 beta(11,3)\n\n\nplt.hist(posterior_sample[\"pi_sample\"], # 绘制 10000次抽样结果\n         edgecolor = \"black\",\n         color=\"grey\",\n         alpha = 0.7,\n         density=True,)\n\nsns.despine()","outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/EAA7BAAB1597442FB63BC0D31EC2D720/s2ld1oz9t7.png\">"],"text/plain":["<Figure size 640x480 with 1 Axes>"]},"metadata":{},"output_type":"display_data"}],"execution_count":7},{"cell_type":"markdown","metadata":{"id":"F52DB418F7524EE2BC8D4EF969F560E6","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"**增加网格数**  \n\n* 这个结果当然过度简化了后验分布，还记得我们刚刚说，只要取的网格越多，对后验的估计会更准确。  \n\n* 那么，现在$\\pi$不取6个值，而是在0~1之间取101个值$\\pi \\in \\{0, 0.01, 0.02, \\ldots, 0.99, 1\\}$  \n\n* 同样的，我们执行上述四个步骤："},{"cell_type":"code","metadata":{"collapsed":false,"id":"E1DC95EA85A1434D8BF884A223944C3A","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"pi_grid = np.linspace(0, 1, 101)            # 生成一个 101 个点，范围在 [0, 1] 之间\n\nprior = st.beta.pdf(pi_grid, 2, 2)          # 生成 Beta(2,2)\n\n\nlikelihood = st.binom.pmf(9, 10, pi_grid)   # 生成二项分布, 参数为 n=10（总试验次数），k=9（成功次数），以及 pi_grid 中的每个概率值\n\n\nposterior = prior * likelihood / np.sum(prior * likelihood)     # 计算后验概率，即先验概率和似然函数的乘积，然后除以归一化常数（分母和）","outputs":[],"execution_count":10},{"cell_type":"code","metadata":{"collapsed":false,"id":"88F3D2F8590545DB8D1AC93F5BBC0239","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"source":"# 画图\nplt.stem(pi_grid, posterior,\n    linefmt='black',\n    bottom=-1)\nplt.ylim(0, 0.05) \nplt.xlabel('pi_grid')\nplt.ylabel('Posterior')\nsns.despine()","outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/B8D45E9FBC1A4648B153C3AD6F05EE00/s2m1ovmznr.png\">"],"text/plain":["<Figure size 640x480 with 1 Axes>"]},"metadata":{},"output_type":"display_data"}],"execution_count":11},{"cell_type":"markdown","metadata":{"id":"103549C5514C4BB6A8E98B9DA62AC59A","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"可以看到，在此次的抽样中，$\\pi$的取值更加多样"},{"cell_type":"code","metadata":{"collapsed":false,"id":"4C2129D0E1F94B3998A03FBC8AB50B29","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"source":"np.random.seed(84735)\n\nposterior_sample = np.random.choice(pi_grid,        # pi_grid 为抽样的候选集合，\n                                    size = 10000,   # 从 posterior 分布中随机抽取 10000 个样本，使用 np.random.choice()\n                                    p=posterior,    # p=posterior 为抽样的概率分布\n                                    replace=True)   # replace=True 表示有放回抽样\n\nposterior_sample = pd.DataFrame({\"pi_sample\": posterior_sample})    # 将抽取的样本存储在 DataFrame 中，列名为 \"pi_sample\"\n\nposterior_sample.value_counts(normalize=True).reset_index()        # 对 posterior_sample 中的样本进行计数，并使用 normalize=True 将计数转换为相对频率\n                                                                   # 然后通过 reset_index() 将结果重新设置为 DataFrame","outputs":[{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th>pi_sample</th>\n","      <th>0</th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>0</th>\n","      <td>0.82</td>\n","      <td>0.0410</td>\n","    </tr>\n","    <tr>\n","      <th>1</th>\n","      <td>0.84</td>\n","      <td>0.0380</td>\n","    </tr>\n","    <tr>\n","      <th>2</th>\n","      <td>0.79</td>\n","      <td>0.0379</td>\n","    </tr>\n","    <tr>\n","      <th>3</th>\n","      <td>0.80</td>\n","      <td>0.0374</td>\n","    </tr>\n","    <tr>\n","      <th>4</th>\n","      <td>0.83</td>\n","      <td>0.0372</td>\n","    </tr>\n","    <tr>\n","      <th>...</th>\n","      <td>...</td>\n","      <td>...</td>\n","    </tr>\n","    <tr>\n","      <th>60</th>\n","      <td>0.37</td>\n","      <td>0.0004</td>\n","    </tr>\n","    <tr>\n","      <th>61</th>\n","      <td>0.35</td>\n","      <td>0.0004</td>\n","    </tr>\n","    <tr>\n","      <th>62</th>\n","      <td>0.41</td>\n","      <td>0.0003</td>\n","    </tr>\n","    <tr>\n","      <th>63</th>\n","      <td>0.39</td>\n","      <td>0.0002</td>\n","    </tr>\n","    <tr>\n","      <th>64</th>\n","      <td>0.38</td>\n","      <td>0.0002</td>\n","    </tr>\n","  </tbody>\n","</table>\n","<p>65 rows × 2 columns</p>\n","</div>"],"text/plain":["    pi_sample       0\n","0        0.82  0.0410\n","1        0.84  0.0380\n","2        0.79  0.0379\n","3        0.80  0.0374\n","4        0.83  0.0372\n","..        ...     ...\n","60       0.37  0.0004\n","61       0.35  0.0004\n","62       0.41  0.0003\n","63       0.39  0.0002\n","64       0.38  0.0002\n","\n","[65 rows x 2 columns]"]},"execution_count":12,"metadata":{},"output_type":"execute_result"}],"execution_count":12},{"cell_type":"markdown","metadata":{"id":"33A9B814C56F4189BDE74EF32327FEB3","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"**抽样结果图示**  \n\n对比一下抽样得到的后验分布图 和 实际的后验分布$Beta(11, 3)$  \n\n相比于上一次，这一次的结果更加接近真实的后验分布"},{"cell_type":"code","metadata":{"collapsed":false,"id":"0FCCCE3701454B8D9A6AF804B855FBB7","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"source":"x_beta = np.linspace(0,1,10000)             # 生成一个 10000 个点，范围在 [0, 1] 之间\ny_beta = st.beta.pdf(x_beta, 11, 3)         # 生成Beta(11,3)\n\nplt.plot(x_beta, y_beta, color='black')     # 绘制Beta(11,3)\n\nplt.hist(posterior_sample[\"pi_sample\"],     # 绘制10000个抽样结果分布\n         edgecolor = \"black\",\n         color=\"grey\",\n         alpha = 0.7,\n         density=True,)\n\nsns.despine()","outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/D21490F1C43341688D0F2D95493C75EB/s2ld1ojz5b.png\">"],"text/plain":["<Figure size 640x480 with 1 Axes>"]},"metadata":{},"output_type":"display_data"}],"execution_count":12},{"cell_type":"markdown","metadata":{"id":"525A3377F4264077B37E759E5EDC42E3","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"## 练习  \n\n> 📃以**Gamma-Poisson模型**为例来练习网格近似  \n\n设$\\lambda$是中国在亚运会期间每天获得金牌的**平均次数**，我们将$\\lambda$的先验分布设为$\\text{Gamma}(3, 1)$  \n\nY为数据，即每天中国实际获得的金牌数量，则：  \n\n$$  \n\\begin{equation}  \n\\begin{split}  \nY_i|\\lambda & \\stackrel{ind}{\\sim} \\text{Pois}(\\lambda) \\\\  \n\\lambda   & \\sim \\text{Gamma}(3, 1)  . \\\\  \n\\end{split}  \n\\tag{6.2}  \n\\end{equation}  \n$$  \n\n如果在第一天中国获得2枚金牌 ($Y_1 = 2$)，在第二天获得金牌的数量为8 ($Y_2 = 8$)  \n\n后验分布:  \n$$  \n\\lambda | ((Y_1,Y_2) = (2,8)) \\sim \\text{Gamma}(13, 3)  \n$$"},{"cell_type":"markdown","metadata":{"id":"A1D9EF1322D3422483696EE7BC4319F3","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"### step1: 选定一系列离散的$\\lambda$值  \n\n1. 在(0,15)这个范围内选择501个$\\lambda$值  \n\n> 注：由$\\lambda$的先验与似然函数分布图我们可以看到，$\\lambda$>15的可能性很小  \n![ImageName](https://www.bayesrulesbook.com/bookdown_files/figure-html/ch6-gam-pois-1.png)"},{"cell_type":"code","metadata":{"collapsed":false,"id":"D616DF16CD664843B6353169056478DC","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"#===========================================================================\n#                            在(0,15)这个范围内选择501个$\\lambda$值\n#                            在...中填入相应的值\n#===========================================================================\nlambda_grid = np.linspace(...,...,...)","outputs":[],"execution_count":13},{"cell_type":"markdown","metadata":{"id":"C8702A2E644E4BED87D5288A3ABCCB70","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"### step2: 在每一个$\\pi$下，计算先验分布，与似然函数的乘积，计算总和并进行归一化"},{"cell_type":"code","metadata":{"collapsed":false,"id":"54660C0B5AD64F5E8D0F50018B72F614","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"# 定义似然函数，由于该似然函数不能通过调用已有的函数获得，因此我们根据得到的公式定义一个\ndef poisson_likelihood(y):\n    lambda_grid = np.linspace(0,15,501)\n    \n    # Calculate the Poisson likelihood for each lambda\n    likelihood_poisson = np.exp(-len(y) * lambda_grid) * np.power(lambda_grid, np.sum(y)) / np.prod([np.math.factorial(val) for val in y])\n\n    return likelihood_poisson","outputs":[],"execution_count":14},{"cell_type":"code","metadata":{"collapsed":false,"id":"A203012C27A74AE5A1760F0C9B807A20","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"#===========================================================================\n#                            请修改 ... 中的值。\n#===========================================================================\nobserved_data = [...,...]                                   # 定义观察到的数据\n\nlikelihood_poisson = poisson_likelihood(...)                # 生成对应的似然值\n\nprior_gamma = st.gamma.pdf(..., a=..., scale=...)           # 定义先验\n\nposterior_lambda = ...                                      # 计算后验概率，即先验概率和似然函数的乘积，然后除以归一化常数（分母和）\n","outputs":[],"execution_count":15},{"cell_type":"markdown","metadata":{"id":"7839E9CA3DE74019AC533DFABE70963C","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"### step4: 得到后验结果($\\lambda$)的分布后，从这个分布中抽样10000次"},{"cell_type":"code","metadata":{"collapsed":false,"id":"94DC994E712745D8AE34CBE681FD7548","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"#===========================================================================\n#                            请修改 ... 中的值。\n#===========================================================================\nnp.random.seed(84735)\nlambda_sample = np.random.choice(...,                   # 输入抽样集合\n                                 size = ...,            # 输入抽样的样本个数\n                                 p=...,                 # 输入抽样的概率分布\n                                 replace=True)          # replace=True 表示有放回抽样\n\nlambda_sample = pd.DataFrame({\"lambda_sample\": lambda_sample})","outputs":[],"execution_count":null},{"cell_type":"code","metadata":{"collapsed":false,"id":"886DED602290414CAC233B53F4881377","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"x_gamma = np.linspace(0,15,10000)                   # 生成一个 10000 个点，范围在 [0, 15] 之间\ny_gamma = st.gamma.pdf(x_gamma, a=13, scale=1/3)    # 生成Gamma(13,3)下, x对应的概率值\n\nplt.plot(x_gamma, y_gamma, color='black')           # 绘制Beta(11,3)\n\nplt.hist(lambda_sample[\"lambda_sample\"],            # 绘制10000个抽样结果分布\n         edgecolor = \"black\",\n         bins=15,\n         color=\"grey\",\n         alpha = 0.7,\n         density=True,)\n\nsns.despine()","outputs":[],"execution_count":16},{"cell_type":"markdown","metadata":{"id":"C9FA2DCF17C0482D99C0A3E9FE440AD5","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"### 网格近似的局限性  \n\n😎或许你也注意到了，上述模型还是只有一个参数的情况  \n\n当模型的参数越来越多时，网格近似就会遇到“维数灾难问题”  \n\n在这里我们只需要简单理解为：当模型参数越来越多时，网格近似需要的计算成本也越来越高。"},{"cell_type":"markdown","metadata":{"id":"D998538A237D4E59A9081DAEDD86AF57","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"## Markov chains via pymc"},{"cell_type":"markdown","metadata":{"id":"B8CCE041668D425B8A349469719818FF","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"**“Markov chain Monte Carlo(MCMC)”**  \n\n在计算具有多个参数的后验分布时，有一种有效的方法叫做马尔科夫链蒙特卡洛(MCMC)，可以用于模拟和**近似**后验概率分布  \n\n- 前一个MC中的Markov 来源于一个俄国数学家的名字 Andrey Markov。  \n  - 他研究并提出一种数学方法，被命名为马尔科夫链。用于描述状态空间中经过从一个状态到另一个状态的转换的随机过程。  \n  - 该过程要求具备 **“无记忆性”**，即下一状态的概率分布只能由当前状态决定，在时间序列中它前面的事件均与之无关。  \n  - 马尔可夫链在机器学习和深度学习中非常广泛，例如分层隐马尔可夫模型HMM (Hidden Markov Model)和卡尔曼滤波模型。  \n- 而后一个MC (Monte Carlo)则是一个赌场的名字。据说来自于 20 世纪 40 年代绝密核武器项目的，其中斯坦尼斯拉夫-乌拉姆、约翰-冯-诺伊曼和他们在洛斯阿拉莫斯国家实验室的合作者使用马尔可夫链来模拟和更好地理解中子旅行现象（Eckhardt，1987）。洛斯阿拉莫斯团队将他们的工作称为 \"蒙特卡洛\"，据说这一选择是受到法国里维埃拉富丽堂皇的蒙特卡洛赌场的启发。  \n\n\n<table>  \n    <tr>  \n        <td><img src=\"https://cdn.kesci.com/upload/s2ld84yupl.jpg?imageView2/0/w/320/h/320\" alt=\"图片1\"></td>  \n        <td><img src=\"https://cdn.kesci.com/upload/s2ld8ix3d7.jpg?imageView2/0/w/320/h/320\" alt=\"图片2\"></td>  \n    </tr>  \n</table>  \n"},{"cell_type":"markdown","metadata":{"id":"5F93FC66F7294B1F92E3693B91723E33","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"\n**MCMC的采样特点**  \n\n1. 和网格法类似，MCMC并不会从后验分布$f(\\theta|y)$中直接采样  \n2. 并且区分与网格法，MCMC的样本不是互相独立的 ———下一个样本值依赖于上一个样本值  \n\n![](https://pic2.zhimg.com/80/v2-696651b37f2dc6a472c314f1ba78194d_1440w.webp)  \n> source: https://zhuanlan.zhihu.com/p/250146007  \n"},{"cell_type":"markdown","metadata":{"id":"1DC2DC7B37C143709116BB9A39A8B612","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"假设后验分布仍是关于$\\theta$的分布，$\\left\\lbrace \\theta^{(1)}, \\theta^{(2)}, \\ldots, \\theta^{(N)} \\right\\rbrace$ 构成了一个长度为N的**马尔科夫链**  \n\n- ——>其中，$\\theta^{(2)}$ 的结果依赖于 $\\theta^{(1)}$，$\\theta^{(3)}$ 的结果依赖于 $\\theta^{(2)}$.... $\\theta^{(3)}$ 的结果依赖于 $\\theta^{(i+1)}$ 的结果依赖于 $\\theta^{(i)}$  \n- ——> 总的来说，$\\theta^{(i+1)}$ 的结果依赖于 $\\theta^{(i)}$和收集到的数据y  \n\n$$  \nf\\left(\\theta^{(i + 1)} \\; | \\; \\theta^{(1)}, \\theta^{(2)}, \\ldots, \\theta^{(i)}, y\\right) = f\\left(\\theta^{(i + 1)} \\; | \\; \\theta^{(i)}, y\\right)  .  \n$$  \n\n> 注意：$\\theta^{(i+1)}$ 依赖于 $\\theta^{(i)}$，而 $\\theta^{(i)}$ 又依赖于$\\theta^{(i-1)}$  \n> 然而，$\\theta^{(i+1)}$ 独立于 $\\theta^{(i-1)}$，这是马尔科夫链 **“无记忆性”**的特点。"},{"cell_type":"markdown","metadata":{"id":"F8F8950EFE0D4EEB9124722A896EC050","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"在lec7中，我们将更详细地介绍MCMC的具体原理，现在有很多的统计工具来帮助我们进行MCMC模拟。  \n\n在这节课中，我们主要介绍如何使用Pymc来进行MCMC模拟的基本过程。"},{"cell_type":"markdown","metadata":{"id":"DD2D9B9C68AA4A1CB6FEBBDEFF01F53F","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"## What is Pymc?  \n\n![Image Name](https://cdn.kesci.com/upload/image/rjvix0jc6b.png?imageView2/0/w/320/h/320)  \n\n\nPyMC是一个在Python中基于**概率编程语音(probabilistic program language, PPL)**的工具，允许用户使用简单的Python API构建贝叶斯模型，并使用马尔可夫链蒙特卡洛（MCMC）方法进行拟合。  \n\n特点：  \n\n- PyMC致力于使贝叶斯建模尽可能简单，使用户能够专注于问题本身，而不是方法。  \n- 包括最先进的推断算法，包括MCMC（NUTS）和变分推断（ADVI）。  \n- 高效快速：使用PyTensor作为计算后端，通过C、Numba或JAX进行编译，将模型运行在GPU上，并受益于复杂的图优化。  \n- 提供后续的分析方法：包括概率分布、高斯过程、ABC、SMC等等。它与ArviZ集成得很好，用于可视化和诊断，以及与Bambi集成用于高级混合效应模型。  \n- 注重社区：在讨论区提问问题，参加MeetUp活动，关注我们的Twitter，并开始贡献。"},{"cell_type":"markdown","metadata":{"id":"8A620C1B2BAE456F925552CFE293ED18","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"## A Beta-Binomial example in pymc  \n\n *我们使用的仍然是本节课一开始的  投票数例子*  \n\n\n> 先验分布为：$\\pi    \\sim \\text{Beta}(2, 2)$  \n> 似然函数为：$Y|\\pi  \\sim \\text{Bin}(10, \\pi)$  \n> 总投票数为10，支持票数为9(Y = 9)  \n\n\n我们尝试使用 Pymc 来表达和设定 Beta-Binomial 模型"},{"cell_type":"code","metadata":{"collapsed":false,"id":"6857A9C520F048828B7E629D2D1B40CF","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"import pymc as pm\nimport arviz as az","outputs":[],"execution_count":17},{"cell_type":"code","metadata":{"collapsed":false,"id":"8FE932B269304C72A7B70E89C130CE4C","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"Y = 9\n\nbb_model = pm.Model()\n\nwith bb_model:\n    pi = pm.Beta('pi', alpha=2, beta=2)\n    likelihood = pm.Binomial('likelihood', n=10, p=pi, observed=Y)","outputs":[],"execution_count":18},{"cell_type":"markdown","metadata":{"id":"07858946BAAE46CD9F110C7013375CAC","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"**代码解释**  \n\n```python  \nY = 9  \n\nbb_model = pm.Model()  \n\nwith bb_model:  \n    pi = pm.Beta('pi', alpha=2, beta=2)  \n    likelihood = pm.Binomial('likelihood', n=10, p=pi, observed=Y)  \n```\n\n在pymc中，一个模型的定义通常包含了先验和似然两部分(其复杂性视模型复杂程度而定)：  \n\n1. 设立容器  \n    * 在pymc中，你需要创建一个`pm.Model()`来容纳你模型中的变量  \n\n    * 接下来你需要定义模型里的各种参数，即`with bb_model:`，表明接下来你对模型中各参数的设定，都会被添加到该模型中  \n\n    * 或者也可以直接一步写成：`with pm.Model() as bb_model:`  \n\n2. 定义先验  \n    * 在这个例子中，我们需要对$\\pi$进行定义  \n\n3. 定义似然  \n    * 我们通过observed = Y，将收集到的数据传入似然函数中  \n"},{"cell_type":"markdown","metadata":{"id":"0D905E1F65AC48DCA6AB1B8739736405","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"### 模型可视化  \n\n可以通过PyMC3自带的可视化工具将模型关系可视化"},{"cell_type":"code","metadata":{"collapsed":false,"id":"1E257042C529452782011CC5259051BF","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"source":"pm.model_to_graphviz(bb_model)","outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/DBA902FA6BDC417091875E940BCDC04F/s2ldvx3gg1.svg\">"],"text/plain":["<graphviz.graphs.Digraph at 0x7fad62a4d1c0>"]},"execution_count":19,"metadata":{},"output_type":"execute_result"}],"execution_count":19},{"cell_type":"markdown","metadata":{"id":"38A92439375D46B18F7ECE7A56CBED55","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"### 使用mcmc进行采样  \n\n在以下例子中，将 MCMC 采样方法得到的参数样本赋值为 `trace`。  \n- 使用 `sample` 方法进行 MCMC 采样模拟过程。  \n- 设置参数 `draws` 来控制 MCMC 采样的次数。  \n- `chains` 表示同时运行几条MCMC链。  \n\n我们可以使用 `arviz` 的方法 `plot_trace` 来可视化该结果  \n- 左图为参数分布图  \n- 右图为 trace 图，代表随着采样的进行(即x轴1-5000次采样)，每个参数值的大小(即y轴为每个采样参数的大小)。"},{"cell_type":"code","metadata":{"collapsed":false,"id":"6FB379634558474995631B2DF775B0D9","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"source":"#采样过程仍在该容器中进行\nwith bb_model:\n    trace = pm.sample(draws=5000,                   # 使用mcmc方法进行采样，draws为采样次数\n                      chains=1,                     # 设置MCMC的链数\n                      random_seed=84735)            # 设置随机状态，以获得与notebook相同的结果\n\naz.plot_trace(trace)","outputs":[{"name":"stderr","output_type":"stream","text":["Auto-assigning NUTS sampler...\n","Initializing NUTS using jitter+adapt_diag...\n","Sequential sampling (1 chains in 1 job)\n","NUTS: [pi]\n"]},{"data":{"text/html":["\n","<style>\n","    /* Turns off some styling */\n","    progress {\n","        /* gets rid of default border in Firefox and Opera. */\n","        border: none;\n","        /* Needs to be in here for Safari polyfill so background images work as expected. */\n","        background-size: auto;\n","    }\n","    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n","        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n","    }\n","    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n","        background: #F44336;\n","    }\n","</style>\n"],"text/plain":["<IPython.core.display.HTML object>"]},"metadata":{},"output_type":"display_data"},{"name":"stderr","output_type":"stream","text":["Sampling 1 chain for 1_000 tune and 5_000 draw iterations (1_000 + 5_000 draws total) took 2 seconds.\n","Only one chain was sampled, this makes it impossible to run some convergence checks\n"]},{"data":{"text/plain":["array([[<Axes: title={'center': 'pi'}>, <Axes: title={'center': 'pi'}>]],\n","      dtype=object)"]},"execution_count":21,"metadata":{},"output_type":"execute_result"},{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/989687C68D51424EB54B92AAD497E105/s2m2xjg75k.png\">"],"text/plain":["<Figure size 1200x200 with 2 Axes>"]},"metadata":{},"output_type":"display_data"}],"execution_count":21},{"cell_type":"markdown","metadata":{"id":"B614F8B6576F46B99D9C0C8328161972","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"### 采样的时间进程  \n\n下图展示了第一条Makov链的前20个采样结果和前200个结果  \n"},{"cell_type":"code","metadata":{"collapsed":false,"id":"AD0F529B87EC41D58F0568598E5080F5","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"source":"# 选取第一条Makov链的前20个采样结果和前200个结果\nsamples_20 = trace.posterior[\"pi\"].sel(chain=0).values[:20]\nsamples_200 = trace.posterior[\"pi\"].sel(chain=0).values[:200]\n\n# 绘图\naz.plot_trace({\"pi\": samples_20})\naz.plot_trace({\"pi\": samples_200})","outputs":[{"data":{"text/plain":["array([[<Axes: title={'center': 'pi'}>, <Axes: title={'center': 'pi'}>]],\n","      dtype=object)"]},"execution_count":22,"metadata":{},"output_type":"execute_result"},{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/E93FBF71046C485E8A1CC81535104D2D/s2m31njrdq.png\">"],"text/plain":["<Figure size 1200x200 with 2 Axes>"]},"metadata":{},"output_type":"display_data"},{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/E93FBF71046C485E8A1CC81535104D2D/s2m31npmnj.png\">"],"text/plain":["<Figure size 1200x200 with 2 Axes>"]},"metadata":{},"output_type":"display_data"}],"execution_count":22},{"cell_type":"code","metadata":{"collapsed":false,"id":"FDD1B73049C143C4AAF0D353F3508D87","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"source":"post_pi = pd.DataFrame({\"pi\": trace.posterior[\"pi\"].values.reshape(-1)})\npost_pi","outputs":[{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th>pi</th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>0</th>\n","      <td>0.741947</td>\n","    </tr>\n","    <tr>\n","      <th>1</th>\n","      <td>0.544469</td>\n","    </tr>\n","    <tr>\n","      <th>2</th>\n","      <td>0.858457</td>\n","    </tr>\n","    <tr>\n","      <th>3</th>\n","      <td>0.415141</td>\n","    </tr>\n","    <tr>\n","      <th>4</th>\n","      <td>0.739664</td>\n","    </tr>\n","    <tr>\n","      <th>...</th>\n","      <td>...</td>\n","    </tr>\n","    <tr>\n","      <th>4995</th>\n","      <td>0.730776</td>\n","    </tr>\n","    <tr>\n","      <th>4996</th>\n","      <td>0.938275</td>\n","    </tr>\n","    <tr>\n","      <th>4997</th>\n","      <td>0.821699</td>\n","    </tr>\n","    <tr>\n","      <th>4998</th>\n","      <td>0.821699</td>\n","    </tr>\n","    <tr>\n","      <th>4999</th>\n","      <td>0.821699</td>\n","    </tr>\n","  </tbody>\n","</table>\n","<p>5000 rows × 1 columns</p>\n","</div>"],"text/plain":["            pi\n","0     0.741947\n","1     0.544469\n","2     0.858457\n","3     0.415141\n","4     0.739664\n","...        ...\n","4995  0.730776\n","4996  0.938275\n","4997  0.821699\n","4998  0.821699\n","4999  0.821699\n","\n","[5000 rows x 1 columns]"]},"execution_count":23,"metadata":{},"output_type":"execute_result"}],"execution_count":23},{"cell_type":"markdown","metadata":{"id":"974EE981E9744CC69460BFF7334D8A2B","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"### 采样结果可视化  \n\n* 将采样结果(5000次采样)对比真实的后验分布(黑线)Beta(11, 3)  \n\n* 可以看到这个采样结果很好地近似了后验分布"},{"cell_type":"code","metadata":{"collapsed":false,"id":"D362F8E853254706A3BF60F949A62835","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"fragment"},"tags":[],"trusted":true},"source":"fig, (axs1, axs2) = plt.subplots(1, 2, figsize=(15, 4))\n\n#绘制采样结果直方图\nsns.histplot(data=post_pi, \n             x=\"pi\", \n             bins=30,\n             ax=axs1,\n             edgecolor='#20699d', \n             color=\"#6497b1\",\n             alpha = 1)\n\n#绘制采样结果密度分布图\nsns.kdeplot(data=post_pi,\n             x=\"pi\",\n             color='#6497b1',\n             fill=True,\n             alpha = 1,\n             ax=axs2)\n\n#绘制真实后验分布图\nx = np.linspace(0.2, 1, 10000)\ny = st.beta.pdf(x, 11, 3)\naxs2.plot(x, y, color='black')\n\n\nsns.despine()","outputs":[{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/2BC4BE8B55604945AF2218F1684C2612/s2m31un5v8.png\">"],"text/plain":["<Figure size 1500x400 with 2 Axes>"]},"metadata":{},"output_type":"display_data"}],"execution_count":24},{"cell_type":"markdown","metadata":{"id":"88EC5879846941698E998279CE88EFDC","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"## 练习  \n> 📃以**Gamma-Poisson模型**为例来练习使用pymc进行MCMC模拟  \n\n设$\\lambda$则是每小时事件发生的平均次数，我们将$\\lambda$的先验分布设为$\\text{Gamma}(3, 1)$  \n\nY是一个小时内事件发生的次数，则：  \n\n$$  \n\\begin{equation}  \n\\begin{split}  \nY_i|\\lambda & \\stackrel{ind}{\\sim} \\text{Pois}(\\lambda) \\\\  \n\\lambda   & \\sim \\text{Gamma}(3, 1)  . \\\\  \n\\end{split}  \n\\tag{6.2}  \n\\end{equation}  \n$$  \n\n如果在第一个小时内事件发生的次数为2 ($Y_1 = 2$)，第二个小时内事件发生的次数为8 ($Y_2 = 8$)  \n\n后验分布:  \n$$  \n\\lambda | ((Y_1,Y_2) = (2,8)) \\sim \\text{Gamma}(13, 3)  \n$$"},{"cell_type":"code","metadata":{"collapsed":false,"id":"D29892D9FCFB40EB970897403DB9B464","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"import pymc as pm\nimport arviz as az","outputs":[],"execution_count":3},{"cell_type":"code","metadata":{"collapsed":false,"id":"512B9E9C88FC4ABD8AD53894D80D291B","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"#===========================================================================\n#                            请修改 ... 中的值。\n#===========================================================================\n\ndata = [...,...]\n#1. 设立容器\n\nwith pm.Model() as gp_model:\n    #2. 定义先验\n    lambda_ = pm.Gamma('lambda_', alpha=..., beta=...)\n    #3. 定义似然\n    obs = pm.Poisson('obs', mu=lambda_, observed=...)","outputs":[],"execution_count":15},{"cell_type":"code","metadata":{"collapsed":false,"id":"5DF28BEE427F4F81B0C04D2ED2DF3676","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"#===========================================================================\n#                            请修改 ... 中的值。\n#===========================================================================\nwith gp_model:\n    trace = pm.sample(draws=...,                    # 例如采样次数设为500-5000次\n                      chains=1,                     # 链为1条\n                      random_seed=84735)","outputs":[],"execution_count":16},{"cell_type":"code","metadata":{"collapsed":false,"id":"40B74FC26203418183EE59E85E43E10C","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"trace","outputs":[],"execution_count":17},{"cell_type":"code","metadata":{"collapsed":false,"id":"AC1BB0B5DEF2457692D4559A81C7E88B","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"az.plot_trace(trace)","outputs":[],"execution_count":21},{"cell_type":"code","metadata":{"collapsed":false,"id":"098AFCD3A3F942089AD001CFBA34F4E1","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"skip"},"tags":[],"trusted":true},"source":"post_lambda = pd.DataFrame({\"lambda\": trace.posterior[\"lambda_\"].values.reshape(-1)})","outputs":[],"execution_count":18},{"cell_type":"code","metadata":{"collapsed":false,"id":"FE166A2DE62142C4AFE9425EE8EF6354","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"fig, (axs1, axs2) = plt.subplots(1, 2, figsize=(15, 4))\n\nsns.histplot(data=post_lambda, \n             x=\"lambda\", \n             bins=30,\n             ax=axs1,\n             edgecolor='#20699d', \n             color=\"#6497b1\",\n             alpha = 1)\n\nsns.kdeplot(data=post_lambda,\n             x=\"lambda\",\n             color='#6497b1',\n             fill=True,\n             alpha = 1,\n             ax=axs2)\n\nx = np.linspace(0, 15, 10000)\ny = st.gamma.pdf(x, a=13, scale=1/3)\naxs2.plot(x, y, color='black')\n\n\nsns.despine()","outputs":[],"execution_count":20},{"cell_type":"markdown","metadata":{"id":"FFDA2C0E1F68410EA32753CE568E1457","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"## 对马尔科夫链的诊断 (Markov chain diagnostics) "},{"cell_type":"markdown","metadata":{"id":"590390D8AAC14506AEDB64A8CB2BE79E","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"在之前的例子中，我们看到马尔科夫链遍历了参数（π 或 λ）的样本空间，并最终模拟了收敛到后验的随机样本。  \n\n而模拟采样方法的关键是 \"近似 \"和 \"收敛\"  \n\n这就引出了以下问题：  \n\n- 好的马尔可夫链是什么样的？  \n- 如何判断马尔可夫链样本是否产生了合理的后验近似值？  \n- 马尔可夫链样本量应该多大？  \n\n诊断的意义：因为我们不知道真实的后验，这意味着我们无法将模拟结果与 \"真实 \"后验结果进行比较。因此，诊断具有重要的参考意义。"},{"cell_type":"markdown","metadata":{"id":"B954DF09676046B78B8262C754503BD1","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"在本节中，我们将重点介绍几种诊断方法  \n\n- 可视化诊断：轨迹图 (trace plot)和平行链 (parallel chains)。  \n- 诊断指标：有效样本大小 ESS(effective sample size)、自相关性 (autocorrelation)和 $\\hat{R}$ (R-hat)。  \n    \n> 这些诊断方法的使用应综合考虑。由于没有一种视觉或数值诊断方法是万能的，因此将它们放在一起考虑时，可以更全面地反映马尔可夫链的质量。"},{"cell_type":"markdown","metadata":{"id":"06C4A0AB059D4CC3868ABEDE9B46A9F6","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"### trace plots  \n\n轨迹图 (trace plot) 是一种可视化方法，用于显示马尔可夫链的轨迹，即马尔可夫链在参数空间中的位置。  \n- 图中的横坐标是时间步长，纵坐标是参数的值。  \n- 轨迹图往往结合参数的后验分布一起展示，以便于观察采样样本对于参数分布的代表性。"},{"cell_type":"markdown","metadata":{"id":"26AC8ADAFB6A41E99072BA7B502DDB3B","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"我们先来看一下正常的轨迹图是什么样：  \n\n\n![Image Name](https://cdn.kesci.com/upload/s2le92xkaq.png?imageView2/0/w/500/h/500)  \n\n- 这是之前 Beta-Binomial model 的示例。  \n- 上图：轨迹图看起来像一堆白噪声，没有明显的趋势或现象，这意味着**链是稳定的**。  \n- 下右：展示了正常的参数后验分布的情况。"},{"cell_type":"markdown","metadata":{"id":"BD1803353EDA445A9621F2E011322B8A","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"一个糟糕的轨迹图可能看起来像这。显然这种条件下的后验参数分布在用于后续推断时是必须非常谨慎的：  \n\n\n![Image Name](https://cdn.kesci.com/upload/s2le9nbz6q.png?imageView2/0/w/450/h/450)  \n\n图中上部分链 A：  \n\n- 链 A 中的轨迹图显示，它在 5000 次迭代后还没有稳定下来，并且它只探索了 0.6 到 0.9 之间的参数值。  \n  - 下降趋势也暗示着链值之间存在着很强的相关性--它们看起来并不像独立的噪音。  \n  - 虽然马尔可夫链本身具有依赖性，但它们越像噪声(独立样本的表现)，得出的后验近似值误差就越小（粗略地说）。  \n- 结合后验分布图：**它的后验近似值高估了真实分布中央部分，而完全低估了这个范围之外值的可信度**。  \n\n图中下部分链 B：   \n- 链 B 表现出不同的问题，迹线图中的部分区域存在**两条完全·平直的线**所显示的那样，  \n  - 这意味着，**当它采样到较小的参数值时，往往会陷入这个值的附近，** 这也表明了一种局部的高相关性。  \n  - 链条 B 在陷入困境时，会对后验参数左侧尾部的值进行**过度采样**。         \n- 结合密度图  \n  - 虽然链 B 后验分布和真实分布的重合性更好，但它存在多峰分布，原因在于过度采样。"},{"cell_type":"markdown","metadata":{"id":"444724E8C0614EBB8DCD9D64290C119B","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"如果我们得到一个糟糕的轨迹图，我们可以采取一些措施：  \n\n- 1.检查模型。确定假定的**先验模型**和**数据模型**是否合适？    \n- 2.对数据链进行**更多迭代**。一些不理想的短期连锁趋势可能会在长期内得到改善。"},{"cell_type":"markdown","metadata":{"id":"9C843B9BAFEB43809A144D2F0AA807FE","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"### 有效样本大小 ESS(effective sample size)和自相关性 (autocorrelation)  \n\nMCMC 链的样本存在相关性，如本身的特性一样，后一个参数的采样过程与前一个参数有关，这导致间隔一个位置的参数样本间的相关性很高。  \n\n然后，MCMC 链又要求“无记忆性”，即样本之间最好是独立的。  \n\n🤔那么如何理解 MCMC链中独立与相关的矛盾？  \n- 需要清楚的是，相关性仅限于间隔一个位置的参数样本。  \n- 而独立性的要求是针对于参数样本间大于一个间隔的情况的。  \n\n自相关就可以很好的描述样本之间的相关性；而有效样本量针对于描述样本的独立性。"},{"cell_type":"markdown","metadata":{"id":"89BD857AE90945AD9F49813AB41D2174","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"**自相关性 (autocorrelation)**  \n\n- 可用于评估马尔科夫链的样之间的相关性。  \n- 强烈的自相关性或依赖性是一件坏事--它**与较小的有效样本比相伴而生**，因此提供了一个警告信号，表明我们得出的后验近似值可能不可靠。  \n    \n根据马尔可夫链的简单构造，链值之间必然存在一些自相关性--一个链值（π(i)）取决于前一个链值（π(i-1)），而前一个链值（π(i-2)）又取决于前一个链值（π(i-3)），依此类推。**这种依赖链还意味着，每个链值都在一定程度上依赖于之前的所有链值**。  \n\n- 例如，π(i) 取决于π(i-1)，而π(i-1) 取决于π(i-2)，因此π(i) 也取决于π(i-2)。  \n- 然而，**这种依赖性或自相关性会逐渐消失**。  \n    \n> 这就像托布勒的地理学第一定律：万事万物都与其他事物相关，**但近处的事物比远处的事物更相关**。"},{"cell_type":"markdown","metadata":{"id":"0F7A4AD85AD240539D75560DD8E42AC4","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"一个自相关良好的例子：  \n\n\n![Image Name](https://cdn.kesci.com/upload/s2lefvdivs.png?imageView2/0/w/640/h/640)  \n\n- 图左是之前提到的轨迹图(trace plot)。  \n- 右图是自相关图，x轴是样本之间的间隔距离，y 轴为自相关性。  \n- 例如 0 代表所有参数采样和自己的相关性 (这显然为1)，1代表所有参数和间隔一个距离参数的样本的相关性，该相关性依然很高，大约为 0.5，表明相差仅 1 步的链值之间存在中等程度的相关性。随后自相关性迅速下降，到间隔 5 个样本时相关性下降为 0。  \n- 也就是说，相差间隔越大，马尔科夫链值之间的相关性越小，这是一个好消息，表明马尔科夫链正在**快速混合**，即在后验可信 π 值范围内快速移动，从而至少模拟了一个独立样本。  \n- 对于间隔解释的例子：如果有100个参数形成的样本，编号为 1-100。 那么间隔10，就代表，10-100组成的样本。为了让 10-100 和原始的样本长度相同，我们可以使用 1-90和 10-100 的样本进行相关分析。需要注意的是，存在更好的数学方法可以避免间隔越大，样本量越小的问题。扩展阅读: https://zhuanlan.zhihu.com/p/77072803  \n  "},{"cell_type":"markdown","metadata":{"id":"10EE71C2B733497D881FC033380D5174","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"一个自相关糟糕的例子：  \n\n\n![Image Name](https://cdn.kesci.com/upload/s2lew3t7qx.png?imageView2/0/w/640/h/640)  \n- 自相关曲线的缓慢下降表明，链值之间的依赖关系不会很快消失。  \n- 相差整整 20 步的马尔可夫链值之间存在大约 0.9 的相关性！  \n- 由于其链值与前值的联系非常紧密，因此该链的**混合速度很慢**，**需要很长时间才能充分探索后验的全部范围**。  \n- 因此，我们应该谨慎使用该链来近似后验。  \n"},{"cell_type":"markdown","metadata":{"id":"F444CCEB0F8F4529AECD1889AFCD67E5","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"**重要概念：Fast vs slow mixing Markov chains**  \n\n- 快速混合：指链表现出与独立样本的行为，即混合链值之间的自相关性迅速下降，有效样本大小比相当大。  \n    \n- 慢速混合：指链不具有独立样本的特征，即混合链值之间的自相关性下降非常缓慢，有效样本大小比很小。换句话说，链在后验可信值范围内 \"缓慢 \"移动。"},{"cell_type":"markdown","metadata":{"id":"C7B4AA9C79314E26840308FF0AD86FCE","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"如果我们的链混合速度太慢该怎么办呢？一些推荐的策略：  \n\n- 第一， \"运行更长的链\"。如果迭代次数足够多，即使是缓慢的混合链最终也能产生良好的后验近似值。  \n    \n- 第二，使得马尔科夫链更薄。  \n    \n    - 例如，在包含 5000 个样本的链中，我们每间隔一个参数来保留样本，并丢弃其余的样本，例如：{π(2), π(4), π(6), . ., π(5000)}.  \n        \n    - 甚至，我们可以每间隔十样本来保留样本：{π(10), π(20), π(30), ..., π(5000)}.  \n        \n- 通过舍弃中间的抽样，我们可以消除低间隔带来的强相关性。  \n    \n    - 还记得吗，π(20) 与稀疏链中的前一个值 (π(10))的相关性小于与原始链中的前一个值 (π(19))的相关性。"},{"cell_type":"markdown","metadata":{"id":"4439ADCB800B41F49D985D04F34E0EEC","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"**一个使 MCMC 链变得更薄(thin)来减少自相关的示例：**  \n\n我们使用 thin=10 来使得链变得比之前薄10倍，也就是从5000个样本中，每隔10个样本宝保留一个，总共500个样本将被保留。  \n\n\n![Image Name](https://cdn.kesci.com/upload/s2lew3t7qx.png?imageView2/0/w/640/h/640)  \n\n![Image Name](https://cdn.kesci.com/upload/s2lf7kf4e8.png?imageView2/0/w/640/h/640)  \n\n- 上图为之前没有变薄的MCMC链。下图为变薄了10倍的链。  \n- 可以看到，自相关性产生了快速的下降。  \n- 但是，这可能任然无法完全消除自相关性。  \n\n🤔一个值得思考的问题：  \n\n- 使用 thin=10 时，我们可用的样本变得更少，是否值得损失 90% 的原始样本值呢？"},{"cell_type":"markdown","metadata":{"id":"1A67573290734829AAF9673973E80191","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"\n**有效样本量比率 Effective sample size ratio**  \n\n为了避免链中相关性的影响，有效样本量 (effective sample size)被提出，以用于描述链中独立样本的数量。因为根据 MCMC 的特性，我们需要链存在无记忆性，也就是大于两个间隔的参数样本间不存在或者存在很小的相关性。有效样本量计算了这些不存在或者存在很小的相关性样本的数量。  \n\n$\\frac{N_{neff}}{N}$  \n\n- N 表示依赖马尔可夫链的实际样本量或长度。  \n- Neff 为马尔可夫链的有效样本大小，量化了产生等效精确后验近似所需的独立样本数量。  \n- Neff 越大越好，但马尔可夫链近似的准确性通常只能达到较小独立样本的准确性。也就是说，  \n- 通常情况下，Neff < N，因此有效样本大小比小于1。我们往往会对有效样本量比小于 0.1 的马尔可夫链产生怀疑，即 Neff/N < 10%。  \n\n假设马尔科夫链的长度为 20,000，并且包括 6800 个独立样本，即有效样本量 ESS = 6800。  \n- 因此，有效样本量比率 EFF=34%，代表相当于我们只使用了 34% 的独立样本。  \n- 由于这个比例高于 0.1，这是是可以接受的。并且 ESS > 400。  \n  \n> 更详细的定义请参见 Vehtari 等人（2021）。如需了解有效样本量, R-hat等指标的的联系，请参阅 Vats 和 Knudson (2018)。"},{"cell_type":"markdown","metadata":{"id":"9B7C931BA10D430CAC4BED5A0690ACE9","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"**补充，如何计算ESS**  \n\n\nESS表示与独立和均匀采样相当的有效样本数量，可以使用以下公式计算：  \n\nESS = N / τ  \n\n- 其中，N是总的采样数量。  \n- τ为自相关时间。  \n\n估计自相关时间（Autocorrelation Time，τ）：自相关时间表示样本之间的相关性消失所需的步数。它可以通过截尾自相关函数进行估计。自相关时间τ可以通过求和截尾自相关函数的面积来估计，可以使用以下公式计算：  \n\nτ = 1 + 2 * Σ(TACF(lag))  \n\n- 其中，lag表示滞后的范围。  \n\n\n计算自相关函数（Autocorrelation Function，ACF）：  \n- 首先，计算采样序列中不同滞后（lag）的自相关系数。  \n- 滞后为0时，自相关系数等于1；滞后为1时，自相关系数表示第一个样本与第二个样本之间的相关性；以此类推。  \n- 自相关函数的计算可以使用统计软件或编程语言中的相关函数。  \n\n计算截尾自相关函数（Truncated Autocorrelation Function，TACF）：  \n- 为了避免计算无限滞后的自相关系数，通常会截取滞后一定范围内的自相关系数。  \n- 选择一个合适的截断范围，例如滞后至自相关系数小于0.1，可以避免过多计算。  \n"},{"cell_type":"markdown","metadata":{"id":"68A28AD449C7418AA21B8DCE97C675CB","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"### R-hat and comparing parallel chains  \n\n自相关和有效样本量仅针对于描述一条MCMC链的有效性。假如我们获得了四条并行的马尔可夫链。我们不仅希望看到每条链的稳定性，还希望看到**四条链的一致性**。  \n\n可以想象，如果四条链的结果并不一致，那么推断结果的可靠性就很难保证。  \n\n同样，我们首先展示一个良好的并行链的结果：  \n\n\n![Image Name](https://cdn.kesci.com/upload/s2lg1juesc.png?imageView2/0/w/640/h/640)  \n\n- 四条链产生的后验近似值几乎没有差别，这证明我们的模拟是稳定的。"},{"cell_type":"markdown","metadata":{"id":"4474CDA0D16A48A49327EC8AEE2A21BC","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"为了便于比较和展示一个糟糕的并行链的结果，我们用更短的马尔可夫链模拟来演示。  \n\n\n![Image Name](https://cdn.kesci.com/upload/s2lfz6esan.png?imageView2/0/w/640/h/640)  \n\n- 运行四条并行链仅进行 100 次迭代，每次产生 50 个样本量：  \n\n虽然这些链的轨迹图表现出相似的随机行为，但它们对应的密度图却不尽相同，因此产生的后验近似值也不尽相同。  \n\n这表明，如果我们只进行了 100 次迭代就停止模拟，那结果是不稳的。"},{"cell_type":"markdown","metadata":{"id":"9BFF894B724B4434A8E7EED0015A4615","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"此外，不同链样本之间的变异性的一致性非常重要。  \n\n之前的结果同样显示，其中四条链具有相似的变异：  \n\n![Image Name](https://cdn.kesci.com/upload/s2lg1juesc.png?imageView2/0/w/640/h/640)  \n\n\n- 中图显示，四条链的可变性几乎是相同的。  \n- 因此，右图中，所有链的变异与综合为单个链的变异性相似。"},{"cell_type":"markdown","metadata":{"id":"DE65C9EBCA9A4AB9915304A31C979D40","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"一个糟糕的例子：  \n\n\n![Image Name](https://cdn.kesci.com/upload/s2lg5bgzbr.png?imageView2/0/w/640/h/640)  \n\n\n- 四条平行链产生了相互冲突的后验近似值（中图）  \n- 当我们将这些链组合在一起时，后验近似值不稳定且较差（右图）。  \n- 因此，所有链中参数值的范围和变异性远大于任何单个链中参数值的范围和变异性。"},{"cell_type":"markdown","metadata":{"id":"7F4BA7CB8692454893FE944A7D87795F","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"我们可以直观地看出，所有链条上的数值变化与单个并行链条内部的数值变化之间的关系非常重要。具体来说  \n\n- 在 \"好的 \"马尔可夫链模拟中，所有并行链的变异性总和将与任何单个链内的变异性**大致相当**；  \n- 在 \"坏的 \"马尔可夫链模拟中，所有并行链的变异性总和可能**超过**每个链内的典型变异性  \n    \n我们可以使用 $\\hat{R}$ 量化综合链变异性和链内变异性之间的关系。"},{"cell_type":"markdown","metadata":{"id":"8BC91CBE67C34B14A263A7F381A69C9E","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"$\\hat{R}$ 指标（称 \"R-hat\"）为并行链的可视化比较提供了具体的数学指标。  \n\n\n$$  \n\\hat{R} = \\frac{V_{combined}}{V_{within}}  \n$$  \n- $\\hat{R}$ 通过比较所有链上采样值的变异性和每个单独链上的变异性来解决一致性问题。  \n- combined 代表所有链的样本值的总变异。  \n- within 代表每个链内部样本值的变异。  \n\n理想情况下，$\\hat{R}$ ≈1，反映了平行链的稳定性。  \n- 相反，$\\hat{R}$ > 1表示链之间的变异不稳定，即组合链内的变异性大于单个链内的变异性。  \n- 虽然不存在黄金法则，但$\\hat{R}$比率大于1.05会给模拟的稳定性带来危险信号。  \n\n> 更详细的定义请参见 Vehtari 等人（2021）。如需了解有效样本量, R-hat等指标的的联系，请参阅 Vats 和 Knudson (2018)。"},{"cell_type":"markdown","metadata":{"id":"5A6B7FAF7143401BAEE0BE25A55C1CED","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"对于示例1  \n\n![Image Name](https://cdn.kesci.com/upload/s2lg1juesc.png?imageView2/0/w/640/h/640)  \n- $\\hat{R}$ 值等于 1，这反映了四个并行链之间和链内部的变异性存在很好的一致性。  \n\n对于示例2  \n\n![Image Name](https://cdn.kesci.com/upload/s2lg5bgzbr.png?imageView2/0/w/640/h/640)  \n\n- $\\hat{R}$ = 5.35，这表明所有链值的方差总和是每个链内典型方差的 5 倍多。这远远超过了 1.05 的标准，充分证明假设平行链并没有产生一致的后验近似值。  \n- 因此 MCMC 的模拟近似是不稳定的。其后验分布的结果不可用于后续推断。"},{"cell_type":"markdown","metadata":{"id":"92203754DBD442018523F463EEEF78A1","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"**代码演示**  \n\n使用 pymc 建模的模型默认使用 MCMC 方法进行参数估计。  \n\n并且我们可以很方便的使用 `arviz` 对MCMC进行可视化和基于数学化标准的诊断"},{"cell_type":"markdown","metadata":{"id":"F6274866224342B1AA33D0C852F97C84","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"首先，我们建立一个简单的 beta 模型，并对该模型使用MCMC采样方法  \n- 设置采样数 draws 为 5000  \n- 设置 MCMC 链的数量为 4，和上面的示例一样。  \n- 设置 MCMC 调试的数量 tune 为0，而默认 tune 为 1000。  \n    - 这里暂时理解 tune 为提前进行 MCMC 采样 1000次，然后丢弃掉这些采样。  \n    - 根据 MCMC 的性质，链中越早的采样越不稳定。因此，丢掉这些采样可以保证后续的采样更加稳定。  \n    - 后续课程会详细介绍该参数的意义。"},{"cell_type":"code","metadata":{"collapsed":false,"id":"58AF7E53055540118E119C5FB3A10385","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"import pymc as pm\nimport arviz as az\n\n# 设置数据\nY = 9\n# 创建PyMC模型\nbb_model = pm.Model()\n\nwith bb_model:\n    pi = pm.Beta('pi', alpha=2, beta=2)\n    likelihood = pm.Binomial('likelihood', n=10, p=pi, observed=Y)\n\n#采样过程仍在该容器中进行\nwith bb_model:\n    trace = pm.sample(draws=5000,                   # 设置MCMC采样次数\n                      tune=0,                       # 设置调试和需要丢弃的采样数\n                      chains=4,                     # 设置4条链\n                      random_seed=84735)","outputs":[{"name":"stderr","output_type":"stream","text":["Auto-assigning NUTS sampler...\n","Initializing NUTS using jitter+adapt_diag...\n","Multiprocess sampling (4 chains in 4 jobs)\n","NUTS: [pi]\n"]},{"data":{"text/html":["\n","<style>\n","    /* Turns off some styling */\n","    progress {\n","        /* gets rid of default border in Firefox and Opera. */\n","        border: none;\n","        /* Needs to be in here for Safari polyfill so background images work as expected. */\n","        background-size: auto;\n","    }\n","    progress:not([value]), progress:not([value])::-webkit-progress-bar {\n","        background: repeating-linear-gradient(45deg, #7e7e7e, #7e7e7e 10px, #5c5c5c 10px, #5c5c5c 20px);\n","    }\n","    .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n","        background: #F44336;\n","    }\n","</style>\n"],"text/plain":["<IPython.core.display.HTML object>"]},"metadata":{},"output_type":"display_data"},{"name":"stderr","output_type":"stream","text":["Sampling 4 chains for 0 tune and 5_000 draw iterations (0 + 20_000 draws total) took 17 seconds.\n"]}],"execution_count":35},{"cell_type":"markdown","metadata":{"id":"44E1CF5F41C64F159C4273551050C1B5","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"绘制 trace 图。  \n\n- 可以看到，每条链的后验分布都重合在一起。  \n- MCMC链的行为表现得类似噪音一样，这是一个好的表现。"},{"cell_type":"code","metadata":{"collapsed":false,"id":"459D1D41D70743F492658F8E8271E8A5","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"az.plot_trace(trace, legend = True)","outputs":[{"data":{"text/plain":["array([[<Axes: title={'center': 'pi'}>, <Axes: title={'center': 'pi'}>]],\n","      dtype=object)"]},"execution_count":39,"metadata":{},"output_type":"execute_result"},{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/AB36F0267D1D43B0AB4BEED33C403D7A/s2m477evz2.png\">"],"text/plain":["<Figure size 1200x200 with 2 Axes>"]},"metadata":{},"output_type":"display_data"}],"execution_count":39},{"cell_type":"markdown","metadata":{"id":"E4C54E5B8FC84CD78D153F2619D5C9EC","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"绘制自相关图  \n\n- 可以看到，随着间隔的下降，自相关也快速下降"},{"cell_type":"code","metadata":{"collapsed":false,"id":"AFFA483BA33149E6B711C5D7FB25563D","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"az.plot_autocorr(trace, max_lag=20, combined=True)","outputs":[{"data":{"text/plain":["<Axes: title={'center': 'pi'}>"]},"execution_count":40,"metadata":{},"output_type":"execute_result"},{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/2C7B7FE8BB1546419986A43DBBFF547F/s2m4axvqog.png\">"],"text/plain":["<Figure size 640x480 with 1 Axes>"]},"metadata":{},"output_type":"display_data"}],"execution_count":40},{"cell_type":"markdown","metadata":{"id":"0174F8881D0E4EE89C0D496B46C40E4D","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"我们使用 arviz 的 `summary`来获得对于 MCMC 链的 R-hat和 有效样本量 ESS 指标。  \n\n- 可以看到 r-hat 接近1，表明链间和各链内的变异的一致性。  \n- 有效样本量 ESS (ess_bulk) 为 4744，有效样本量比率 4744/5000 = 0.95，远大于 0.1。"},{"cell_type":"code","metadata":{"collapsed":false,"id":"CBE1BD73F54C418683F5446BAF0B7C9B","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"az.summary(trace, kind = \"diagnostics\")","outputs":[{"data":{"text/html":["<div>\n","<style scoped>\n","    .dataframe tbody tr th:only-of-type {\n","        vertical-align: middle;\n","    }\n","\n","    .dataframe tbody tr th {\n","        vertical-align: top;\n","    }\n","\n","    .dataframe thead th {\n","        text-align: right;\n","    }\n","</style>\n","<table border=\"1\" class=\"dataframe\">\n","  <thead>\n","    <tr style=\"text-align: right;\">\n","      <th></th>\n","      <th>mcse_mean</th>\n","      <th>mcse_sd</th>\n","      <th>ess_bulk</th>\n","      <th>ess_tail</th>\n","      <th>r_hat</th>\n","    </tr>\n","  </thead>\n","  <tbody>\n","    <tr>\n","      <th>pi</th>\n","      <td>0.001</td>\n","      <td>0.001</td>\n","      <td>4744.0</td>\n","      <td>4366.0</td>\n","      <td>1.0</td>\n","    </tr>\n","  </tbody>\n","</table>\n","</div>"],"text/plain":["    mcse_mean  mcse_sd  ess_bulk  ess_tail  r_hat\n","pi      0.001    0.001    4744.0    4366.0    1.0"]},"execution_count":42,"metadata":{},"output_type":"execute_result"}],"execution_count":42},{"cell_type":"markdown","metadata":{"id":"2EC4C6552F6D4E31A7B9A362C5779966","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"虽然示例的 trace 不存在强烈的自相关，但我们演示如何通过 thin 参数来使得采样更薄，从而减少链的自相关。  \n\n- 我们设置 thin 为10，即每隔10保留一个参数，这导致样本数量只有之前的 1/10，即500  \n- 可以看到，自相关几乎不存在了。"},{"cell_type":"code","metadata":{"collapsed":false,"id":"FB1646A24DFA4CDA9D25B51773A00889","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"thinned_trace = trace.sel(draw=slice(None, None, 10))\n\n# 注意，thinned_trace 只有 500个采样\nthinned_trace.posterior","outputs":[{"data":{"text/html":["<div><svg style=\"position: absolute; width: 0; height: 0; overflow: hidden\">\n","<defs>\n","<symbol id=\"icon-database\" viewBox=\"0 0 32 32\">\n","<path d=\"M16 0c-8.837 0-16 2.239-16 5v4c0 2.761 7.163 5 16 5s16-2.239 16-5v-4c0-2.761-7.163-5-16-5z\"></path>\n","<path d=\"M16 17c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n","<path d=\"M16 26c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n","</symbol>\n","<symbol id=\"icon-file-text2\" viewBox=\"0 0 32 32\">\n","<path d=\"M28.681 7.159c-0.694-0.947-1.662-2.053-2.724-3.116s-2.169-2.030-3.116-2.724c-1.612-1.182-2.393-1.319-2.841-1.319h-15.5c-1.378 0-2.5 1.121-2.5 2.5v27c0 1.378 1.122 2.5 2.5 2.5h23c1.378 0 2.5-1.122 2.5-2.5v-19.5c0-0.448-0.137-1.23-1.319-2.841zM24.543 5.457c0.959 0.959 1.712 1.825 2.268 2.543h-4.811v-4.811c0.718 0.556 1.584 1.309 2.543 2.268zM28 29.5c0 0.271-0.229 0.5-0.5 0.5h-23c-0.271 0-0.5-0.229-0.5-0.5v-27c0-0.271 0.229-0.5 0.5-0.5 0 0 15.499-0 15.5 0v7c0 0.552 0.448 1 1 1h7v19.5z\"></path>\n","<path d=\"M23 26h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n","<path d=\"M23 22h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n","<path d=\"M23 18h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n","</symbol>\n","</defs>\n","</svg>\n","<style>/* CSS stylesheet for displaying xarray objects in jupyterlab.\n"," *\n"," */\n","\n",":root {\n","  --xr-font-color0: var(--jp-content-font-color0, rgba(0, 0, 0, 1));\n","  --xr-font-color2: var(--jp-content-font-color2, rgba(0, 0, 0, 0.54));\n","  --xr-font-color3: var(--jp-content-font-color3, rgba(0, 0, 0, 0.38));\n","  --xr-border-color: var(--jp-border-color2, #e0e0e0);\n","  --xr-disabled-color: var(--jp-layout-color3, #bdbdbd);\n","  --xr-background-color: var(--jp-layout-color0, white);\n","  --xr-background-color-row-even: var(--jp-layout-color1, white);\n","  --xr-background-color-row-odd: var(--jp-layout-color2, #eeeeee);\n","}\n","\n","html[theme=dark],\n","body[data-theme=dark],\n","body.vscode-dark {\n","  --xr-font-color0: rgba(255, 255, 255, 1);\n","  --xr-font-color2: rgba(255, 255, 255, 0.54);\n","  --xr-font-color3: rgba(255, 255, 255, 0.38);\n","  --xr-border-color: #1F1F1F;\n","  --xr-disabled-color: #515151;\n","  --xr-background-color: #111111;\n","  --xr-background-color-row-even: #111111;\n","  --xr-background-color-row-odd: #313131;\n","}\n","\n",".xr-wrap {\n","  display: block !important;\n","  min-width: 300px;\n","  max-width: 700px;\n","}\n","\n",".xr-text-repr-fallback {\n","  /* fallback to plain text repr when CSS is not injected (untrusted notebook) */\n","  display: none;\n","}\n","\n",".xr-header {\n","  padding-top: 6px;\n","  padding-bottom: 6px;\n","  margin-bottom: 4px;\n","  border-bottom: solid 1px var(--xr-border-color);\n","}\n","\n",".xr-header > div,\n",".xr-header > ul {\n","  display: inline;\n","  margin-top: 0;\n","  margin-bottom: 0;\n","}\n","\n",".xr-obj-type,\n",".xr-array-name {\n","  margin-left: 2px;\n","  margin-right: 10px;\n","}\n","\n",".xr-obj-type {\n","  color: var(--xr-font-color2);\n","}\n","\n",".xr-sections {\n","  padding-left: 0 !important;\n","  display: grid;\n","  grid-template-columns: 150px auto auto 1fr 20px 20px;\n","}\n","\n",".xr-section-item {\n","  display: contents;\n","}\n","\n",".xr-section-item input {\n","  display: none;\n","}\n","\n",".xr-section-item input + label {\n","  color: var(--xr-disabled-color);\n","}\n","\n",".xr-section-item input:enabled + label {\n","  cursor: pointer;\n","  color: var(--xr-font-color2);\n","}\n","\n",".xr-section-item input:enabled + label:hover {\n","  color: var(--xr-font-color0);\n","}\n","\n",".xr-section-summary {\n","  grid-column: 1;\n","  color: var(--xr-font-color2);\n","  font-weight: 500;\n","}\n","\n",".xr-section-summary > span {\n","  display: inline-block;\n","  padding-left: 0.5em;\n","}\n","\n",".xr-section-summary-in:disabled + label {\n","  color: var(--xr-font-color2);\n","}\n","\n",".xr-section-summary-in + label:before {\n","  display: inline-block;\n","  content: '►';\n","  font-size: 11px;\n","  width: 15px;\n","  text-align: center;\n","}\n","\n",".xr-section-summary-in:disabled + label:before {\n","  color: var(--xr-disabled-color);\n","}\n","\n",".xr-section-summary-in:checked + label:before {\n","  content: '▼';\n","}\n","\n",".xr-section-summary-in:checked + label > span {\n","  display: none;\n","}\n","\n",".xr-section-summary,\n",".xr-section-inline-details {\n","  padding-top: 4px;\n","  padding-bottom: 4px;\n","}\n","\n",".xr-section-inline-details {\n","  grid-column: 2 / -1;\n","}\n","\n",".xr-section-details {\n","  display: none;\n","  grid-column: 1 / -1;\n","  margin-bottom: 5px;\n","}\n","\n",".xr-section-summary-in:checked ~ .xr-section-details {\n","  display: contents;\n","}\n","\n",".xr-array-wrap {\n","  grid-column: 1 / -1;\n","  display: grid;\n","  grid-template-columns: 20px auto;\n","}\n","\n",".xr-array-wrap > label {\n","  grid-column: 1;\n","  vertical-align: top;\n","}\n","\n",".xr-preview {\n","  color: var(--xr-font-color3);\n","}\n","\n",".xr-array-preview,\n",".xr-array-data {\n","  padding: 0 5px !important;\n","  grid-column: 2;\n","}\n","\n",".xr-array-data,\n",".xr-array-in:checked ~ .xr-array-preview {\n","  display: none;\n","}\n","\n",".xr-array-in:checked ~ .xr-array-data,\n",".xr-array-preview {\n","  display: inline-block;\n","}\n","\n",".xr-dim-list {\n","  display: inline-block !important;\n","  list-style: none;\n","  padding: 0 !important;\n","  margin: 0;\n","}\n","\n",".xr-dim-list li {\n","  display: inline-block;\n","  padding: 0;\n","  margin: 0;\n","}\n","\n",".xr-dim-list:before {\n","  content: '(';\n","}\n","\n",".xr-dim-list:after {\n","  content: ')';\n","}\n","\n",".xr-dim-list li:not(:last-child):after {\n","  content: ',';\n","  padding-right: 5px;\n","}\n","\n",".xr-has-index {\n","  font-weight: bold;\n","}\n","\n",".xr-var-list,\n",".xr-var-item {\n","  display: contents;\n","}\n","\n",".xr-var-item > div,\n",".xr-var-item label,\n",".xr-var-item > .xr-var-name span {\n","  background-color: var(--xr-background-color-row-even);\n","  margin-bottom: 0;\n","}\n","\n",".xr-var-item > .xr-var-name:hover span {\n","  padding-right: 5px;\n","}\n","\n",".xr-var-list > li:nth-child(odd) > div,\n",".xr-var-list > li:nth-child(odd) > label,\n",".xr-var-list > li:nth-child(odd) > .xr-var-name span {\n","  background-color: var(--xr-background-color-row-odd);\n","}\n","\n",".xr-var-name {\n","  grid-column: 1;\n","}\n","\n",".xr-var-dims {\n","  grid-column: 2;\n","}\n","\n",".xr-var-dtype {\n","  grid-column: 3;\n","  text-align: right;\n","  color: var(--xr-font-color2);\n","}\n","\n",".xr-var-preview {\n","  grid-column: 4;\n","}\n","\n",".xr-index-preview {\n","  grid-column: 2 / 5;\n","  color: var(--xr-font-color2);\n","}\n","\n",".xr-var-name,\n",".xr-var-dims,\n",".xr-var-dtype,\n",".xr-preview,\n",".xr-attrs dt {\n","  white-space: nowrap;\n","  overflow: hidden;\n","  text-overflow: ellipsis;\n","  padding-right: 10px;\n","}\n","\n",".xr-var-name:hover,\n",".xr-var-dims:hover,\n",".xr-var-dtype:hover,\n",".xr-attrs dt:hover {\n","  overflow: visible;\n","  width: auto;\n","  z-index: 1;\n","}\n","\n",".xr-var-attrs,\n",".xr-var-data,\n",".xr-index-data {\n","  display: none;\n","  background-color: var(--xr-background-color) !important;\n","  padding-bottom: 5px !important;\n","}\n","\n",".xr-var-attrs-in:checked ~ .xr-var-attrs,\n",".xr-var-data-in:checked ~ .xr-var-data,\n",".xr-index-data-in:checked ~ .xr-index-data {\n","  display: block;\n","}\n","\n",".xr-var-data > table {\n","  float: right;\n","}\n","\n",".xr-var-name span,\n",".xr-var-data,\n",".xr-index-name div,\n",".xr-index-data,\n",".xr-attrs {\n","  padding-left: 25px !important;\n","}\n","\n",".xr-attrs,\n",".xr-var-attrs,\n",".xr-var-data,\n",".xr-index-data {\n","  grid-column: 1 / -1;\n","}\n","\n","dl.xr-attrs {\n","  padding: 0;\n","  margin: 0;\n","  display: grid;\n","  grid-template-columns: 125px auto;\n","}\n","\n",".xr-attrs dt,\n",".xr-attrs dd {\n","  padding: 0;\n","  margin: 0;\n","  float: left;\n","  padding-right: 10px;\n","  width: auto;\n","}\n","\n",".xr-attrs dt {\n","  font-weight: normal;\n","  grid-column: 1;\n","}\n","\n",".xr-attrs dt:hover span {\n","  display: inline-block;\n","  background: var(--xr-background-color);\n","  padding-right: 10px;\n","}\n","\n",".xr-attrs dd {\n","  grid-column: 2;\n","  white-space: pre-wrap;\n","  word-break: break-all;\n","}\n","\n",".xr-icon-database,\n",".xr-icon-file-text2,\n",".xr-no-icon {\n","  display: inline-block;\n","  vertical-align: middle;\n","  width: 1em;\n","  height: 1.5em !important;\n","  stroke-width: 0;\n","  stroke: currentColor;\n","  fill: currentColor;\n","}\n","</style><pre class='xr-text-repr-fallback'>&lt;xarray.Dataset&gt;\n","Dimensions:  (chain: 4, draw: 500)\n","Coordinates:\n","  * chain    (chain) int64 0 1 2 3\n","  * draw     (draw) int64 0 10 20 30 40 50 60 ... 4940 4950 4960 4970 4980 4990\n","Data variables:\n","    pi       (chain, draw) float64 0.6946 0.9428 0.8367 ... 0.8636 0.8551 0.8676\n","Attributes:\n","    created_at:                 2023-10-16T08:13:14.811553\n","    arviz_version:              0.14.0\n","    inference_library:          pymc\n","    inference_library_version:  5.6.0\n","    sampling_time:              17.21974229812622\n","    tuning_steps:               0</pre><div class='xr-wrap' style='display:none'><div class='xr-header'><div class='xr-obj-type'>xarray.Dataset</div></div><ul class='xr-sections'><li class='xr-section-item'><input id='section-be902888-5ff3-47e1-93c3-501ebef71f42' class='xr-section-summary-in' type='checkbox' disabled ><label for='section-be902888-5ff3-47e1-93c3-501ebef71f42' class='xr-section-summary'  title='Expand/collapse section'>Dimensions:</label><div class='xr-section-inline-details'><ul class='xr-dim-list'><li><span class='xr-has-index'>chain</span>: 4</li><li><span class='xr-has-index'>draw</span>: 500</li></ul></div><div class='xr-section-details'></div></li><li class='xr-section-item'><input id='section-9430c1d8-034d-423f-972d-d9d34f5ddcb5' class='xr-section-summary-in' type='checkbox'  checked><label for='section-9430c1d8-034d-423f-972d-d9d34f5ddcb5' class='xr-section-summary' >Coordinates: <span>(2)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>chain</span></div><div class='xr-var-dims'>(chain)</div><div class='xr-var-dtype'>int64</div><div class='xr-var-preview xr-preview'>0 1 2 3</div><input id='attrs-8d9bdac9-7ae7-45c7-ae2e-486ec3a8e8ce' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-8d9bdac9-7ae7-45c7-ae2e-486ec3a8e8ce' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-1ad83eb9-c0fb-49d0-acbf-8064403ca6ac' class='xr-var-data-in' type='checkbox'><label for='data-1ad83eb9-c0fb-49d0-acbf-8064403ca6ac' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([0, 1, 2, 3])</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>draw</span></div><div class='xr-var-dims'>(draw)</div><div class='xr-var-dtype'>int64</div><div class='xr-var-preview xr-preview'>0 10 20 30 ... 4960 4970 4980 4990</div><input id='attrs-27bbf80b-7a92-4163-be7c-bc1cd156d8ee' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-27bbf80b-7a92-4163-be7c-bc1cd156d8ee' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-883b8692-3d84-41b5-b0a8-68787fa7f71d' class='xr-var-data-in' type='checkbox'><label for='data-883b8692-3d84-41b5-b0a8-68787fa7f71d' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([   0,   10,   20, ..., 4970, 4980, 4990])</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-d46de329-da21-4ce8-aa85-d944cf8da38a' class='xr-section-summary-in' type='checkbox'  checked><label for='section-d46de329-da21-4ce8-aa85-d944cf8da38a' class='xr-section-summary' >Data variables: <span>(1)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span>pi</span></div><div class='xr-var-dims'>(chain, draw)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>0.6946 0.9428 ... 0.8551 0.8676</div><input id='attrs-8dae5a22-282a-4971-a354-d1a08b1b5b83' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-8dae5a22-282a-4971-a354-d1a08b1b5b83' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-3745f322-b08d-4f0c-8279-ec46ed91d2cf' class='xr-var-data-in' type='checkbox'><label for='data-3745f322-b08d-4f0c-8279-ec46ed91d2cf' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([[0.69459658, 0.94278139, 0.83667045, ..., 0.72543572, 0.77218586,\n","        0.79339494],\n","       [0.83451846, 0.93718211, 0.79875923, ..., 0.73197361, 0.74135104,\n","        0.67311654],\n","       [0.91612668, 0.87065615, 0.90718409, ..., 0.53466908, 0.88057703,\n","        0.90598199],\n","       [0.81972277, 0.77083496, 0.76215542, ..., 0.86357788, 0.85514079,\n","        0.86764909]])</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-0495ad1f-4f5e-45e2-82c6-f2d5a042905e' class='xr-section-summary-in' type='checkbox'  ><label for='section-0495ad1f-4f5e-45e2-82c6-f2d5a042905e' class='xr-section-summary' >Indexes: <span>(2)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-index-name'><div>chain</div></div><div class='xr-index-preview'>PandasIndex</div><div></div><input id='index-f9904e6d-e7da-460b-8c14-7ed5cfe546b3' class='xr-index-data-in' type='checkbox'/><label for='index-f9904e6d-e7da-460b-8c14-7ed5cfe546b3' title='Show/Hide index repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-index-data'><pre>PandasIndex(Int64Index([0, 1, 2, 3], dtype=&#x27;int64&#x27;, name=&#x27;chain&#x27;))</pre></div></li><li class='xr-var-item'><div class='xr-index-name'><div>draw</div></div><div class='xr-index-preview'>PandasIndex</div><div></div><input id='index-52c2f504-ec7c-4a9d-b8be-97f3b775ad0f' class='xr-index-data-in' type='checkbox'/><label for='index-52c2f504-ec7c-4a9d-b8be-97f3b775ad0f' title='Show/Hide index repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-index-data'><pre>PandasIndex(Int64Index([   0,   10,   20,   30,   40,   50,   60,   70,   80,   90,\n","            ...\n","            4900, 4910, 4920, 4930, 4940, 4950, 4960, 4970, 4980, 4990],\n","           dtype=&#x27;int64&#x27;, name=&#x27;draw&#x27;, length=500))</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-ca29e9a4-1164-4fdb-ad1f-f09cbff34baf' class='xr-section-summary-in' type='checkbox'  checked><label for='section-ca29e9a4-1164-4fdb-ad1f-f09cbff34baf' class='xr-section-summary' >Attributes: <span>(6)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><dl class='xr-attrs'><dt><span>created_at :</span></dt><dd>2023-10-16T08:13:14.811553</dd><dt><span>arviz_version :</span></dt><dd>0.14.0</dd><dt><span>inference_library :</span></dt><dd>pymc</dd><dt><span>inference_library_version :</span></dt><dd>5.6.0</dd><dt><span>sampling_time :</span></dt><dd>17.21974229812622</dd><dt><span>tuning_steps :</span></dt><dd>0</dd></dl></div></li></ul></div></div>"],"text/plain":["<xarray.Dataset>\n","Dimensions:  (chain: 4, draw: 500)\n","Coordinates:\n","  * chain    (chain) int64 0 1 2 3\n","  * draw     (draw) int64 0 10 20 30 40 50 60 ... 4940 4950 4960 4970 4980 4990\n","Data variables:\n","    pi       (chain, draw) float64 0.6946 0.9428 0.8367 ... 0.8636 0.8551 0.8676\n","Attributes:\n","    created_at:                 2023-10-16T08:13:14.811553\n","    arviz_version:              0.14.0\n","    inference_library:          pymc\n","    inference_library_version:  5.6.0\n","    sampling_time:              17.21974229812622\n","    tuning_steps:               0"]},"execution_count":59,"metadata":{},"output_type":"execute_result"}],"execution_count":59},{"cell_type":"code","metadata":{"collapsed":false,"id":"37D10DFC14534CFDB09E5060CC42DEDC","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[],"trusted":true},"source":"az.plot_autocorr(trace, max_lag=10, combined=True)\naz.plot_autocorr(thinned_trace, max_lag=10, combined=True)","outputs":[{"data":{"text/plain":["<Axes: title={'center': 'pi'}>"]},"execution_count":55,"metadata":{},"output_type":"execute_result"},{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/3A8622959CAB4146AEB6D8DEC76F4A6F/s2m56fnvhw.png\">"],"text/plain":["<Figure size 640x480 with 1 Axes>"]},"metadata":{},"output_type":"display_data"},{"data":{"text/html":["<img src=\"https://cdn.kesci.com/upload/rt/3A8622959CAB4146AEB6D8DEC76F4A6F/s2m56fe2t5.png\">"],"text/plain":["<Figure size 640x480 with 1 Axes>"]},"metadata":{},"output_type":"display_data"}],"execution_count":55},{"cell_type":"markdown","metadata":{"id":"7CE63A83CCDC48B694A365B9E85D3652","jupyter":{},"notebookId":"6536255793c31faf0a5a8dc8","runtime":{"status":"default","execution_status":null,"is_visible":false},"scrolled":false,"slideshow":{"slide_type":"slide"},"tags":[]},"source":"### 总结  \n\n随着我们的贝叶斯模型越来越复杂，比如，一个回归模型加入越来越多的预测变量 (预测纣王的当选的可能)  \n\n这些模型的**后验将变得难以直接获得，或者无法获得**。  \n\n我们学习了在这种情况下逼近后验的两种模拟技术：  \n- 网格逼近和马尔科夫链蒙特卡罗。  \n- 两种技术都会产生 N 个 θ 参数值样本。注意，这里的 θ 可能包含多个参数。  \n\n最后，我们学到了一些用于检查 MCMC 模拟质量的诊断方法。  \n- 可视化：直观地检查多个并行链的轨迹图和密度图，以确定模拟的稳定性和混合性。  \n- 数学指标：用有效样本量比、自相关性和 $\\hat{R}$等数值诊断来补充这些视觉诊断。"}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python","nbconvert_exporter":"python","file_extension":".py","version":"3.5.2","pygments_lexer":"ipython3"}},"nbformat":4,"nbformat_minor":2}